1 Introduction

Researchers have previous suggested that the response diversity of a community be measured by the diversity of responses to environmental change. For example, one can measure the response of each of the species’ intrinsic growth rate to temperature, quantify the strength and direction of these responses (e.g., as the first derivative of the response curve), and calculate the diversity of responses (e.g., by calculating variation in the first derivatives among the species in a community). When responses are nonlinear, the response diversity will be a function of the environmental state (i.e. the first derivative is a function of the value of the environmental state). So far we demonstrated this approach for quantifying response diversity in the context of a single environmental factor, but given that multiple environmental factors can change simultaneously, we need an approach that works in that context.

2 The principle

Owen learned about the mathematical principles from these youtube videos:

Imagine that the growth rate of a population depends on two environmental factors, e.g. temperature and salinity. We can represent the dependency as \(G = f(T, S)\), where \(G\) is growth rate, \(T\) is temperature, and \(S\) is salinity. It may be that the dependencies are linear, nonlinear, and with an interaction between temperature and salinity, hence our approach needs to be able to accommodate this phenomena.

The response of growth rate to change in temperature and salinity is the gradient / slope of this surface, with units of growth rate [per time] per temperature [degrees C] per salinity [parts per thousand]. Because the slope (first derivative) of the surface can (when dependencies are nonlinear) vary across the surface (location on the surface), and can vary in different directions on the surface, to calculate a slope we must specify the current environment (location on the surface) and the direction of change in the environment. The location on the curve is the current environmental condition, \((T_0, S_0)\), and the direction of environmental change is the unit vector \(\hat{u} = \langle U_T, U_S \rangle\).

Put another way, we calculate a directional derivative at a point on the response surface. We can write this as \(D_{\hat{u}}f(T_0, S_0)\) and can calculate it as \(f_T(T_0, S_0)U_T + f_S(T_0, S_0)U_S\), where \(f_T\) is the partial derivative of \(f(T, S)\) with respect to \(T\) and \(f_S\) is the partial derivative of \(f(T, S)\) with respect to \(S\).

Efficient evaluating in \(n\) dimensions can be done by taking the dot product of the partial derivatives at the location and the direction unit vector: \(D_{\hat{u}}f(T_0, S_0) = \triangledown f \cdot \hat{u}\) where, \(\triangledown f = \langle f_T, f_S \rangle\). (In R, the dot product of a and b is sum(a*b))

Figure 2.1 is an illustration of the principle of directional derivatives on a surface.

This figure needs considerable improvement! It is currently a keynote illustration.

Figure 2.1: This figure needs considerable improvement! It is currently a keynote illustration.

3 A simulated empirical example

Numerous mathematical functions have been used to represent how organismal performance changes with an environmental driver citation require. Moreover, multiple mathematical functions have been used to represent an interactive effect of two or more environmental drivers on species performance e.g. Thomas et al 2017. We have to decide if we are going to make simulations with different of these functions. This might increase our confidence that the method for calculating response diversity is robust to variation in types of response curve. If we decide not to, then we should likey argue that we think its robust.

3.1 Simulating performance curves

Let us use the Eppley performance curve, which was used, for example, in this paper Bernhardt et al. 2018.

With one environmental variable, the performance (i.e., rate) is given by:

  • \(rate(E) = ae^{bE}(1 - (\frac{E - z}{w/2})^2)\)
  • \(E\) is the values of the environmental factor.
  • \(z\) controls location of maximum.
  • \(w\) controls range of \(E\) over which the rate is positive.
  • \(a\) scaling constant.
  • \(b\) controls rate of increase towards the maximum rate, as \(E\) increases.

Adding a second environmental variable gives:

\(rate(E_1, E_2) = a_1e^{b_1E_1}(1 - (\frac{E_1 - z_1}{w_1/2})^2) + a_2e^{b_2E_2}(1 - (\frac{E_2 - z_2}{w_2/2})^2)\)

In this case, it is clear the effect of \(E_1\) and \(E_2\) is defined as being additive. For example, the value of \(E_2\) does not affect the value of \(E_1\) at which the rate is maximised (\(z_1\)), and vice-versa (see also Figure 3.1)

Nonlinear and additive dependence of a rate on two environmental variables. (a) The value of \(E_1\) at which the rate is maximised is independent of the value of \(E_2\). (b) The value of \(E_2\) at which the rate is maximised is independent of the value of \(E_1\).

Figure 3.1: Nonlinear and additive dependence of a rate on two environmental variables. (a) The value of \(E_1\) at which the rate is maximised is independent of the value of \(E_2\). (b) The value of \(E_2\) at which the rate is maximised is independent of the value of \(E_1\).

Including an interaction. One way to do this is to make the value of \(E_1\) at which the rate is maximised depend on the value of \(E_2\):

\(rate(E_1, E_2) = a_1e^{b_1E_1}(1 - (\frac{(E_1 + z_{int21}*E_2- z_1)}{w_1/2})^2 + a_2e^{b_2E_2}(1 - (\frac{E_2 - z_2}{w_2/2})^2\)

When \(z_{int21} = 0\) then this equation becomes the previously mentioned additive one. When \(z_{int} \neq 0\) then the value of \(E_1\) at which the rate is maximised is a function of the value of \(E_2\). We used this method for adding an interaction due to its simplicity. Other methods could be used, and if also or otherwise used could add confidence about the robustness of the method for calculating response diversity.

Nonlinear and additive dependence of a rate on two environmental variables. (a) The value of \(E_1\) at which the rate is maximised is independent of the value of \(E_2\). (b) The value of \(E_2\) at which the rate is maximised is independent of the value of \(E_1\).

Figure 3.2: Nonlinear and additive dependence of a rate on two environmental variables. (a) The value of \(E_1\) at which the rate is maximised is independent of the value of \(E_2\). (b) The value of \(E_2\) at which the rate is maximised is independent of the value of \(E_1\).

3.2 Simulating multiple species’ performance curves

3.2.1 No interacting environmental effects

First we create (or import) a table of parameter values of each species, with species in the rows and parameters in the columns. In the following example, only values of the \(z\) parameters differ among the species (which determine the location of the maximum rate).

For convenience we then convert the table of parameters into a list-column. We can then easily make performance curves of each of the species, and put those into a list-column in the same table.

## convert parameter table to a list-column of a tibble
par_list <- Partable_2_parlist(par_table)
## add performance curves
species_pars <- tibble(species = paste0("s", 1:s), pars = par_list) %>%
  rowwise() %>%
  mutate(expt = Make_expt(E1_series, E2_series, pars))

Here are some examples of the species’ performance curves (only with additive effects of \(E_1\) and \(E_2\)).

Performance curves for a species with maximum growth at low values of \(E_1\). Without interacting environmental effects.

Figure 3.3: Performance curves for a species with maximum growth at low values of \(E_1\). Without interacting environmental effects.

Performance curves for a species with maximum growth at high values of \(E_1\). Without interacting environmental effects.

Figure 3.4: Performance curves for a species with maximum growth at high values of \(E_1\). Without interacting environmental effects.

Performance curves for a species with maximum growth at low values of \(E_2\). Without interacting environmental effects.

Figure 3.5: Performance curves for a species with maximum growth at low values of \(E_2\). Without interacting environmental effects.

Performance curves for a species with maximum growth at high values of \(E_2\). Without interacting environmental effects.

Figure 3.6: Performance curves for a species with maximum growth at high values of \(E_2\). Without interacting environmental effects.

3.2.2 Interacting environmental effects

And now with interacting environmental effects…

For convenience we then convert the table of parameters into a list-column. We can then easily make performance curves of each of the species, and put those into a list-column in the same table.

## convert parameter table to a list-column of a tibble
par_list <- Partable_2_parlist(par_table)
## add performance curves
species_pars <- tibble(species = paste0("s", 1:s), pars = par_list) %>%
  rowwise() %>%
  mutate(expt = Make_expt(E1_series, E2_series, pars))

Here are some examples of the species’ performance curves (with interacting effects of \(E_1\) and \(E_2\)).

Performance curves for a species with maximum growth at low values of \(E_1\). With interacting environmental effects.

Figure 3.7: Performance curves for a species with maximum growth at low values of \(E_1\). With interacting environmental effects.

Performance curves for a species with maximum growth at high values of \(E_1\). With interacting environmental effects.

Figure 3.8: Performance curves for a species with maximum growth at high values of \(E_1\). With interacting environmental effects.

Performance curves for a species with maximum growth at low values of \(E_2\). With interacting environmental effects.

Figure 3.9: Performance curves for a species with maximum growth at low values of \(E_2\). With interacting environmental effects.

Performance curves for a species with maximum growth at high values of \(E_2\). With interacting environmental effects.

Figure 3.10: Performance curves for a species with maximum growth at high values of \(E_2\). With interacting environmental effects.

3.3 Fitting GAMs to noisy rate observations

Try with and without an interaction. Therefore make two species, one with no interaction z_int = 0 and the other with z_int = 0.1. All other parameters are the same. Note that noise is added to the rate observations.

Bottom line is that the gam picks up an interaction when we have included one in the parameters used to generate the rates, and does not pick one up when we have not. This confirms that our more mechanistic thinking and methods are matching our statistical thinking and methods, and confirms that each are promising, so far.

GAM notes:

  • Owen watched this youtube video Introduction to Generalized Additive Models with R and mgcv by Gavin Simpson (author of the gratia package that we’ve been using to calculate derivatives of GAMs.). By the way, there is a part of this video about Model Selection, about Confidence Intervals, and about p-values for smooths (about 1h:57m to 2h:26) which I think is less useful for our purposes.
  • In the video mentioned above at about 54m:30s, there is the recommendation to estimate the penalisation parameter (\(\lambda\)) via restricted maximum likelihood, which is done by giving the the argument method = "REML in the gam() call (the default in mgcv version 1.8-36 is method="GCV.Cp").
  • Default smoother is a low rank thin plate spline (bs = "tp")
  • `s(E1, E2) will fit a bivariate isotropic thin plate spline. Isoptropic means there is a single smoothness parameter for the smooth. It is sensitive to the scale of \(E1\) and \(E2\).
  • Tensor products (te(E1, E2) have separate marginal basis, and therefore separate smoothness parameters. They are invariant to scales of \(E1\) and \(E2\). Therefore we should be using tensor products.
  • Specifying a model with te(E1, E2) does not then allow inspection of the main effects and interaction term–there is essentially only one smoother which contains the main effects and interaction. So it is difficult to know if the interaction term is required. An alternative is to specify the model as s(E1) + s(E2) + ti(E1, E2), in which the third term is then only the interaction part and allows a decomposition into the different effects. Useful for examining if the interaction is required.
  • To use cubic splines use bs = "cr".
  • There are many smoothers in mgcv
  • Choosing k (the basis complexity, which is the maximum wiggliness) is a bit of an art. The penalty then shrinks the used basis complexity, reducing the effective degrees of freedom. Must check that large enough k is provided; use gam.check(). Only cost to large \(k\) is computational effort.

3.3.1 Without interaction

Performance curves for a species without interacting environmental effects and with some noise in the rate.

Figure 3.11: Performance curves for a species without interacting environmental effects and with some noise in the rate.

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## rate ~ ti(E1) + ti(E2) + ti(E1, E2)
## 
## Parametric coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.0154285  0.0004991   30.91   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##             edf Ref.df         F p-value    
## ti(E1)    3.999  4.000 14145.143  <2e-16 ***
## ti(E2)    3.937  3.997   609.233  <2e-16 ***
## ti(E1,E2) 1.001  1.002     0.023   0.883    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.958   Deviance explained = 95.8%
## -REML = -5815.8  Scale est. = 0.000648  n = 2601

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## rate ~ s(E1) + s(E2)
## 
## Parametric coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.0154285  0.0003952   39.04   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##         edf Ref.df       F p-value    
## s(E1) 8.970  9.000 10199.8  <2e-16 ***
## s(E2) 7.146  8.187   475.4  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.974   Deviance explained = 97.4%
## -REML = -6407.3  Scale est. = 0.00040624  n = 2601

3.3.2 With interaction

Performance curves for a species with interacting environmental effects and with some noise in the rate.

Figure 3.12: Performance curves for a species with interacting environmental effects and with some noise in the rate.

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## rate ~ ti(E1) + ti(E2) + ti(E1, E2)
## 
## Parametric coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.0656082  0.0004674   140.4   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##             edf Ref.df      F p-value    
## ti(E1)    3.999  4.000 5449.0  <2e-16 ***
## ti(E2)    3.954  3.999  523.7  <2e-16 ***
## ti(E1,E2) 7.252  8.924  878.7  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.924   Deviance explained = 92.5%
## -REML =  -5976  Scale est. = 0.00056817  n = 2601

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## rate ~ s(E1) + s(E2)
## 
## Parametric coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.0656082  0.0009042   72.56   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##         edf Ref.df     F p-value    
## s(E1) 8.561  8.945 672.0  <2e-16 ***
## s(E2) 4.421  5.434 101.2  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.717   Deviance explained = 71.8%
## GCV = 0.0021379  Scale est. = 0.0021264  n = 2601

4 Partial derivatives - no interaction

First step in calculating directional derivatives is estimating the two partial derivatives \(f_{E1}(E1_0, E2_0)\) and \(f_{E2}(E1_0, E2_0)\) (please review the section The principle if necessary).

4.1 Getting the partial derivatives

Partial derivatives. Draw response surface for sp 1 and calculate partial derivatives at a specific location (E1 = 300, E2 = 20). To calculate the partial derivative with respect to E1, E2 must be held constant.

Response surface of sp1. The two solid lines show at which level of E1 and E2 each partial derivative is going to be calculated (E1 = 300, E2 = 20).

Figure 4.1: Response surface of sp1. The two solid lines show at which level of E1 and E2 each partial derivative is going to be calculated (E1 = 300, E2 = 20).

Visualising the partial effect of E1 at a fixed level of E2.

Partial effect of E1 on the growth rate of sp 1 when E2 is held constant at E2 = 20.

Figure 4.2: Partial effect of E1 on the growth rate of sp 1 when E2 is held constant at E2 = 20.

Partial derivative with respect to E1 when E2 is constant at 20.

Partial derivative with respect to E1 when E2 is constant at 20.

Figure 4.3: Partial derivative with respect to E1 when E2 is constant at 20.

Partial derivatives with respect to E2 (E1 held constant)

Response surface of sp1. The two solid lines show at which level of E1 and E2 each partial derivative is going to be calculated.

Figure 4.4: Response surface of sp1. The two solid lines show at which level of E1 and E2 each partial derivative is going to be calculated.

Partial effect of E2 on the growth rate of sp 1 when E1 is held constat at E2 = 300

Partial effect of E2 on the growth rate of sp 1 when E1 is held constant at E1 = 300

Figure 4.5: Partial effect of E2 on the growth rate of sp 1 when E1 is held constant at E1 = 300

Partial derivative with respect to E2 when E1 is constant at 300

Partial derivative with respect to E1 when E2 is constant at 20.

Figure 4.6: Partial derivative with respect to E1 when E2 is constant at 20.

Plot the two partial derivatives and relative effects

Summary plot sp1. (a) response surface of sp 1. (b) Partial effect of E1 on the growth rate of sp 1 when E2 is held constant at E2 = 20. (c) Partial derivative with respect to E1 when E2 is constant at 20. (d) Partial effect of E2 on the growth rate of sp 1 when E1 is held constant at E1 = 300. (e) Partial derivative with respect to E2 when E1 is constant at 300.

Figure 4.7: Summary plot sp1. (a) response surface of sp 1. (b) Partial effect of E1 on the growth rate of sp 1 when E2 is held constant at E2 = 20. (c) Partial derivative with respect to E1 when E2 is constant at 20. (d) Partial effect of E2 on the growth rate of sp 1 when E1 is held constant at E1 = 300. (e) Partial derivative with respect to E2 when E1 is constant at 300.

5 Directional derivatives

5.1 No direction of environmental change specified

5.1.1 One point

We start showing how directional derivatives can be calculated even when the direction of the environmental change is unknown. This may be the case when we want to calculate response diversity for future scenarios, and the future direction of environmental change is thus not known. Or we may have data for a species or a community at only one environmental location (E1 = x, E2 = y). It is therefore important to be able to measure directional derivatives when the direction of the environmental change is unknown, as this can provide useful information on response diversity nonetheless, for instance, by taking the mean of the slopes calculates in all directions.

Here, we calculate, for a specific point (E1 = 300, E2 = 20), directional derivatives in all directions.

Directional derivatives calculated in all possible direction for a specific point on the response surface of sp1. Clearly, the slope of the directional derivative depends on the direction (red positive, blue negative). Note: the size of the radius was only chosen for representation purposes, and does not have any implication. The slope of the segments departing from the point have each their fixed slopes independently of the size of the radius.

Figure 5.1: Directional derivatives calculated in all possible direction for a specific point on the response surface of sp1. Clearly, the slope of the directional derivative depends on the direction (red positive, blue negative). Note: the size of the radius was only chosen for representation purposes, and does not have any implication. The slope of the segments departing from the point have each their fixed slopes independently of the size of the radius.

5.1.2 Several points

We can measure all possible directional derivatives also for several points on the surface. This might be the case when we know that a species or a community occurs at multiple locations on the surface (multiple environmental conditions), but we do not know the direction of change.

Directional derivatives calculated in all possible direction for several points on the response surface of sp1. Clearly, the slope of the directional derivative depends on the direction (red positive, blue negative).

Figure 5.2: Directional derivatives calculated in all possible direction for several points on the response surface of sp1. Clearly, the slope of the directional derivative depends on the direction (red positive, blue negative).

5.1.3 Grid of points

Finally, we might do the same for a grid of points on the surface. We may want to do that when we do not have information on where a species or a community is living within the surface, but we know the range of values of E1 and E2.

Directional derivatives calculated in all possible direction for a grid of points on the response surface of sp1. Clearly, the slope of the directional derivative depends on the direction (red positive, blue negative).

Figure 5.3: Directional derivatives calculated in all possible direction for a grid of points on the response surface of sp1. Clearly, the slope of the directional derivative depends on the direction (red positive, blue negative).

6 Partial derivatives - interactive environmental effect

Repeating the same steps to see if it works with interactive environmental effect.

6.1 Getting the partial derivatives

Partial derivative with respect to E1 (E2 is held constant)

Partial derivatives. Draw response surface for sp 2 and calculate partial derivatives at a specific location (E1 = 300, E2 = 20). To calculate the partial derivative with respect to E1, E2 must be held constant.

Response surface of sp1. The two solid lines show at which level of E1 and E2 each partial derivative is going to be calculated.

Figure 6.1: Response surface of sp1. The two solid lines show at which level of E1 and E2 each partial derivative is going to be calculated.

Partial effect of E1 on the growth rate of sp 2 when E2 is held constant at E2 = 20.

Figure 6.2: Partial effect of E1 on the growth rate of sp 2 when E2 is held constant at E2 = 20.

Partial derivative with respect to E1 when E2 is constant at 20

Partial derivative with respect to E1 when E2 is constant at 20.

Figure 6.3: Partial derivative with respect to E1 when E2 is constant at 20.

Partial derivative with respect to E2 (E1 is held constant)

Response surface of sp2. The two solid lines show at which level of E1 and E2 each partial derivative is going to be calculated.

Figure 6.4: Response surface of sp2. The two solid lines show at which level of E1 and E2 each partial derivative is going to be calculated.

Partial effect of E2 on the growth rate of sp 2 when E1 is held constat at E2 = 300.

Partial effect of E2 on the growth rate of sp 2 when E1 is held constant at E1 = 300

Figure 6.5: Partial effect of E2 on the growth rate of sp 2 when E1 is held constant at E1 = 300

Partial derivative with respect to E2 when E1 is constant at 300.

Partial derivative with respect to E1 when E2 is constant at 20.

Figure 6.6: Partial derivative with respect to E1 when E2 is constant at 20.

Plot the two partial derivatives and relative effects

Summary plot sp1. (a) response surface of sp 1. (b) Partial effect of E1 on the growth rate of sp 1 when E2 is held constant at E2 = 20. (c) Partial derivative with respect to E1 when E2 is constant at 20. (d) Partial effect of E2 on the growth rate of sp 1 when E1 is held constant at E1 = 300. (e) Partial derivative with respect to E2 when E1 is constant at 300.

Figure 6.7: Summary plot sp1. (a) response surface of sp 1. (b) Partial effect of E1 on the growth rate of sp 1 when E2 is held constant at E2 = 20. (c) Partial derivative with respect to E1 when E2 is constant at 20. (d) Partial effect of E2 on the growth rate of sp 1 when E1 is held constant at E1 = 300. (e) Partial derivative with respect to E2 when E1 is constant at 300.

7 Directional derivatives

7.1 No direction of environmental change specified

We start showing how directional derivatives can be calculated even when the direction of the environmental change is unknown. Again, this may be the case when we only know one environmental location in which the species or community is present.

In this case, we calculate, for a specific point (E1 = 300, E2 = 20), directional derivatives in all directions.

7.1.1 One point

Directional derivatives calculated in all possible direction for a specific point on the response surface of sp2. Clearly, the slope of the directional derivative depends on the direction (red positive, blue negative).

Figure 7.1: Directional derivatives calculated in all possible direction for a specific point on the response surface of sp2. Clearly, the slope of the directional derivative depends on the direction (red positive, blue negative).

7.1.2 Several points

We can measure all possible directional derivatives also for several points on the surface. This might be the case when we have multiple locations on the surface (multiple environmental conditions), but we do not know the direction of change.

Directional derivatives calculated in all possible direction for several points on the response surface of sp1. Clearly, the slope of the directional derivative depends on the direction (red positive, blue negative).

Figure 7.2: Directional derivatives calculated in all possible direction for several points on the response surface of sp1. Clearly, the slope of the directional derivative depends on the direction (red positive, blue negative).

7.1.3 Grid of points

Finally, we might do the same for a grid of points on the surface.

Directional derivatives calculated in all possible direction for a grid of points on the response surface of sp1. Clearly, the slope of the directional derivative depends on the direction (red positive, blue negative).

Figure 7.3: Directional derivatives calculated in all possible direction for a grid of points on the response surface of sp1. Clearly, the slope of the directional derivative depends on the direction (red positive, blue negative).

8 Response diversity calculation

Environmental variables may show different correlations between each other. The increase in one environmental variable may be directly correlated with the increase of another one (positive correlation), or vice versa, the increase in one driver may be correlated to a decrease in the other one (negative correlation). Yet, two environmental variable may change over time, or space, completely independently.
We may imagine that these different types of relationships between two environmental variables could determine specific trends in response diversity.

To explore this hypothesis, we calculate now response diversity for two communities composed of 3 spp in 4 different cases: 1. Unknown direction of the environmental change 2. Direction of env change is given by the time series,and E1 and E2 change over time independently 3. Direction of env change is given by the time series,and E1 and E2 change over time with positive correlation 4. Direction of env change is given by the time series,and E1 and E2 change over time with negative correlation

We want to see if any consistent trend appears in the two communities when E1 and E2 have different correlations.

[Only additives effects of the drivers on sp responses - we may change this later, so using an example with interactive effect]

Steps:

  1. Simulate spp performance curves with the modified Eppley function.

  2. Fit response surface for each sp (done with GAMs)

  1. Data wrangling and partials derivatives calculations

8.0.1 Unknown direction of the environmental change

Table showing the calculated response diversity for one of the two communities when the direction of the environmental change is unknown. In this case, directional derivatives for each point where calculated in all directions and then the mean was used to calculate response diversity (only first 6 rows shown).

E1_ref E2_ref s4 s6 s11 rdiv sign Med
273.15 0 0e+00 9e-07 4e-07 1.000000 0.0945711 1.000001
273.15 1 -1e-07 6e-07 4e-07 1.000000 0.1681298 1.000001
273.15 2 -1e-07 2e-07 5e-07 1.000000 0.2156433 1.000001
273.15 3 -1e-07 -2e-07 5e-07 1.000000 0.5279619 1.000001
273.15 4 -1e-07 -5e-07 5e-07 1.000001 0.9810117 1.000001
273.15 5 -1e-07 -9e-07 5e-07 1.000001 0.7519762 1.000001

[Plotting response diversity - think about how to plot RD when we do not know the direction of env change. Surface?]

8.0.2 E1 and E2 change independently over time

This example mimics a situation where the two environmental variables change over time completely independently. This is a common situation in field studies, where multiple drivers of environmental change are not correlated one another.

In this case the direction of the environmental change is given by the change of E1 and E2 over time.

Time series of E1 and E2 changing independently over time.

Figure 8.1: Time series of E1 and E2 changing independently over time.

Table showing the calculated response diversity for one of the two communities when the two environmental variables change independently over time (only first 6 rows shown).

time E1_ref E2_ref s11 s4 s6 rdiv sign Med
1 298.1471 25.01125 -1.00e-06 0.0e+00 -3.60e-06 1.000002 0.0170768 1.000009
2 298.1470 25.01726 1.47e-05 -1.8e-06 1.06e-05 1.000007 0.2174546 1.000009
3 298.1428 25.02150 -3.40e-06 5.0e-07 -3.00e-07 1.000002 0.2503919 1.000009
4 298.1507 24.98623 3.00e-07 -1.0e-07 -2.80e-06 1.000001 0.1665656 1.000009
5 298.1494 25.01227 -2.40e-06 4.0e-07 6.00e-07 1.000001 0.3891664 1.000009
6 298.1510 25.00318 1.82e-05 -2.1e-06 1.73e-05 1.000009 0.2058180 1.000009

Plot response diversity over time

(ref:RDindependentplot) Directional derivatives and response diversity with known direction of env change. E1 and E2 change independently over time. a and b: Species directional derivatives over time. c and d: Response diversity measured as similarity-based diversity metric. e and : Response diversity measured as divergence (sign sensitive).

(ref:RDindependentplot)

Figure 8.2: (ref:RDindependentplot)

8.0.3 E1 and E2 change with negative correlation

This example mimics a situation where the two environmental variables change over time with negative correlation. This is common in field studies, where one environmental variable (e.g. CO2 concentration in oceans) increases, while another (e.g. pH) decreases e.g. Shirayama & Thornton (2005).

Creating a time series with E1 and E2 changing over time with negative correlation.

Time series of E1 and E2 changing with negative correlation over time.

Figure 8.3: Time series of E1 and E2 changing with negative correlation over time.

Table showing the calculated response diversity for one of the two communities when the two environmental variables change with negative correlation over time (only first 6 rows shown).

time E1_ref E2_ref s11 s4 s6 rdiv sign Med
1 297.8961 25.25392 1.17e-05 -1.4e-06 9.50e-06 1.000006 0.2199884 1.000009
2 297.2093 25.94070 2.90e-06 -4.0e-07 5.70e-06 1.000003 0.1396403 1.000009
3 296.2748 26.87522 8.50e-06 -9.0e-07 -2.00e-07 1.000004 0.1915512 1.000009
4 298.8003 24.34969 2.31e-05 -2.8e-06 1.40e-05 1.000012 0.2136118 1.000009
5 298.2508 24.89917 1.64e-05 -2.0e-06 1.14e-05 1.000008 0.2163878 1.000009
6 295.3336 27.81635 1.91e-05 -2.1e-06 5.40e-06 1.000009 0.2014375 1.000009

Plot response diversity over time for the two communities

(ref:RDnegativeplot) Directional derivatives and response diversity with known direction of env change. E1 and E2 change with negative correlation over time. a and b: Species directional derivatives over time. c and d: Response diversity measured as similarity-based diversity metric. e and : Response diversity measured as divergence (sign sensitive).

(ref:RDnegativeplot)

Figure 8.4: (ref:RDnegativeplot)

8.0.4 E1 and E2 change with positive correlation

Finally, two environmental variables can show positive correlation over time. A typical example is given by the positive correlation between air temperature and UV radiation e.g. Häder at al. 2015.

Let us create a time series with E1 and E2 changing over time with positive correlation

Time series of E1 and E2 changing with positive correlation over time.

Figure 8.5: Time series of E1 and E2 changing with positive correlation over time.

Table showing the calculated response diversity for one of the two communities when the two environmental variables change with positive correlation over time (only first 6 rows shown).

time E1_ref E2_ref s11 s4 s6 rdiv sign Med
1 297.9946 23.86346 2.43e-05 -2.9e-06 4.70e-06 1.000012 0.2142890 1.00001
2 297.7153 25.07388 1.22e-05 -1.4e-06 1.34e-05 1.000007 0.1835128 1.00001
3 296.5400 22.72664 2.09e-05 -2.5e-06 -1.54e-05 1.000016 0.8498835 1.00001
4 297.8392 24.50819 4.50e-06 -6.0e-07 -6.50e-06 1.000005 0.8110268 1.00001
5 297.9687 25.20803 -2.06e-05 2.4e-06 -1.77e-05 1.000010 0.2060611 1.00001
6 300.8728 27.62332 4.44e-05 -5.1e-06 9.70e-06 1.000022 0.2069973 1.00001

Plot response diversity over time

(ref:RDpositiveplot) Directional derivatives and response diversity with known direction of env change for community 1 and 2. E1 and E2 change with negative correlation over time. a and b: Species directional derivatives over time. c and d: Response diversity measured as similarity-based diversity metric. e and : Response diversity measured as divergence (sign sensitive).

(ref:RDpositiveplot)

Figure 8.6: (ref:RDpositiveplot)

Now, we visualize the relationship between different correlations between the two environmental variables and response diversity.

(ref:plotcorrelations) Correlation types and response diversity. a and c: correlation types and response diversity measured as dissimilarity in the first derivatives (sign insensitive) for community 1 and 2 respectively. c and d. correlation types and response diversity measured as divergence in the first derivatives (sign sensitive) for community 1 and 2 respectively

(ref:plotcorrelations)

Figure 8.7: (ref:plotcorrelations)

It is hard to derive any conclusion from this plot…

9 Empirical example

We use data coming from an experiment where individual ciliates species have been exposed to a gradient of nutrient, light, and their combinations in a factorial design. We first show how to calculate the partial derivatives, then we calculate the directional derivatives based on a simulated time series (in the original experiment, the level of the treatments have been kept constant throughout the expt duration). Finally, we assemble random composed communities and calculate response diversity for each of them.

9.0.1 Load data set and look at species responses

(ref:ciliates) Species responses to the environmental drivers. a. Species responses to nutrient concentrations. b. Species responses to light intensity.

(ref:ciliates)

Figure 9.1: (ref:ciliates)

9.0.2 Fittig GAMs on empirical data

species E1 E2 predicted
Loxocephallus 0.55 1.00 1255.552
Loxocephallus 0.55 1.02 1248.094
Loxocephallus 0.55 1.04 1240.636
Loxocephallus 0.55 1.06 1233.178
Loxocephallus 0.55 1.08 1225.720
Loxocephallus 0.55 1.10 1218.260

9.0.3 Plotting surface for a sp

Figure 9.2: Response surface fitted with GAM. High non-linearity. Need to understand how to fix axis’ values!

10 Partial derivatives for a single species

10.1 E1 - Nutrients

First, we calculate the partial derivative with respect to nutrient concentration keeping light intensity constant at 5.

Response surface of Colpidium. The two solid lines show at which level of nutrients and light each partial derivative is going to be calculated. Not sure we get the gray areas…

Figure 10.1: Response surface of Colpidium. The two solid lines show at which level of nutrients and light each partial derivative is going to be calculated. Not sure we get the gray areas…

Partial effect of nutrient concentration on the growth rate of Colpidium when light intensity is held constant at 5.

Figure 10.2: Partial effect of nutrient concentration on the growth rate of Colpidium when light intensity is held constant at 5.

Partial derivative with respect to nutrient when light intensity is constant at 5.

Figure 10.3: Partial derivative with respect to nutrient when light intensity is constant at 5.

10.2 E2 - Light intensity

Second, we calculate the partial derivative with respect to light intensity keeping nutrient concentration constant at 2.67.

Response surface of Colpidium. The two solid lines show at which level of nutrients and light each partial derivative is going to be calculated. Not sure we get the gray areas…

Figure 10.4: Response surface of Colpidium. The two solid lines show at which level of nutrients and light each partial derivative is going to be calculated. Not sure we get the gray areas…

Partial effect of light intensity on the growth rate of Colpidium when nutrient concentration is held constant at 2.67.

Figure 10.5: Partial effect of light intensity on the growth rate of Colpidium when nutrient concentration is held constant at 2.67.

Partial derivative with respect to nutrient when light intensity is constant at 5.

Figure 10.6: Partial derivative with respect to nutrient when light intensity is constant at 5.

10.2.1 Plot surface and partial derivatives

Plot the two partial derivatives and relative effects

Summary plot of Colpidium (a) response surface of Colpidium (b) Partial effect of nutrient concentration on the density of Colpidium when light intensity is held constant at 5. (c) Partial derivative with respect to nutrient concentration when light intensity is held constant at 5. (d) Partial effect of light intensity on the growth rate of Colpidium when nutrient concentration is held constant at 2.67. (e) Partial derivative with respect to light intensity when nutrient concentration is held constant at 2.67.

Figure 10.7: Summary plot of Colpidium (a) response surface of Colpidium (b) Partial effect of nutrient concentration on the density of Colpidium when light intensity is held constant at 5. (c) Partial derivative with respect to nutrient concentration when light intensity is held constant at 5. (d) Partial effect of light intensity on the growth rate of Colpidium when nutrient concentration is held constant at 2.67. (e) Partial derivative with respect to light intensity when nutrient concentration is held constant at 2.67.

11 Directional deriviatives

To calculate the directional derivatives for all spp used in the experiment, we first create a time series with nutrient concentration and light intensity changing randomly over time, we fit GAMs individually for each species, and then we calculate partial derivatives.

Time series of nutrient concentration and light intensity changing over time.

Time series of (a) nutrient concentration and (b) light intensity changing over time.

Figure 11.1: Time series of (a) nutrient concentration and (b) light intensity changing over time.

Table with calculated partial derivatives for each sp at different times (only first 6 rowa shown).

sp time E1_ref E2_ref pd_E1 pd_E2
Loxocephallus 1 4.41 2.34 -1400.4493 -1292.1978
Loxocephallus 2 1.15 3.44 -1230.7261 -1036.6754
Loxocephallus 3 2.91 3.08 -1681.4269 -984.1988
Loxocephallus 4 0.61 1.48 188.3544 2584.2893
Loxocephallus 5 4.65 4.04 -511.4889 382.0343
Loxocephallus 6 4.85 3.94 -622.4786 580.1252

11.1 Calculating response diversity for a specific community composition

First, we need to calculate the directional derivatives in the direction of the env change.

sp time E1_ref E2_ref pd_E1 pd_E2 nxt_value_E1 nxt_value_E2 del_E1 del_E2 unit_vec_mag uv_E1 uv_E2 dir_deriv
Loxocephallus 1 4.41 2.34 -1400.4493 -1292.1978 1.15 3.44 -3.26 1.10 3.4405813 -0.9475143 0.3197134 913.8127
Loxocephallus 2 1.15 3.44 -1230.7261 -1036.6754 2.91 3.08 1.76 -0.36 1.7964409 0.9797149 -0.2003962 -998.0149
Loxocephallus 3 2.91 3.08 -1681.4269 -984.1988 0.61 1.48 -2.30 -1.60 2.8017851 -0.8209052 -0.5710645 1942.3331
Loxocephallus 4 0.61 1.48 188.3544 2584.2893 4.65 4.04 4.04 2.56 4.7828025 0.8446930 0.5352510 1542.3451
Loxocephallus 5 4.65 4.04 -511.4889 382.0343 4.85 3.94 0.20 -0.10 0.2236068 0.8944272 -0.4472136 -628.3405
Loxocephallus 6 4.85 3.94 -622.4786 580.1252 1.07 6.88 -3.78 2.94 4.7887368 -0.7893522 0.6139406 847.5173

Then we can calculate response diversity for an ipotetical community contaning all the species tested i this experiment.

time E1_ref E2_ref Coleps_irchel Colpidium Euplotes Paramecium_bursaria Paramecium_caudatum Stylonychia2 rdiv sign Med
1 4.41 2.34 21.631135 289.21778 0.0766887 50.72726 6.993002 14.47996 5.995332 0.0000000 5.986313
2 1.15 3.44 2.118059 60.47016 -0.1122998 -36.86137 -118.478575 -24.40440 5.805878 0.6758378 5.986313
3 2.91 3.08 25.917109 202.91030 -0.1339960 185.73265 77.681554 -13.88633 5.999998 0.1281047 5.986313
4 0.61 1.48 -18.270294 87.35641 -0.1367671 107.91723 -92.974506 -20.24055 5.755278 0.9256180 5.986313
5 4.65 4.04 -44.464615 -395.88035 0.1787494 -93.59877 -235.794177 -38.59871 5.994347 0.0009026 5.986313
6 4.85 3.94 41.038765 403.86375 -0.1372904 119.07608 237.985367 34.07151 5.998117 0.0006797 5.986313

Plot response diversity over time

(ref:RDempiricalplot) Directional derivatives and response diversity with known direction of env change. a. Species directional derivatives over time. b. Response diversity measured as similarity-based diversity metric. c. Response diversity measured as divergence (sign sensitive).

(ref:RDempiricalplot)

Figure 11.2: (ref:RDempiricalplot)

11.2 Different community compositions

Now we calculate response diversity for different community compositions and we comapre them.

Plotting

(ref:RD-compositions) Directional derivatives and response diversity with known direction of env change for three different communities. a. Species directional derivatives over time. b. Response diversity measured as similarity-based diversity metric. c. Response diversity measured as divergence (sign sensitive).

(ref:RD-compositions)

Figure 11.3: (ref:RD-compositions)

11.3 To do

  • The mathematical formulation of the interaction term is not very good. Too much changes when the \(z_int\) term is changed.
  • Collect data on response of growth rate to temperature and salinity for a collection of species. [We will simulate this, including response surfaces with combinations of linear and nonlinear responses, and without and with interactions among temperature and salinity.]
  • Fit a response surface to that data. [Do with GAMs with interactions.]
  • Specify the species present in a community and the temperature and salinity environment it occurs in. [We will simulate this, with different directions of change in temperature and salinity, e.g. independent, positively correlated, negatively correlated.]
  • Show examples of response diversity for the simulated environmental change scenarios and community compositions. [ Make function to calculate directional derivative at location, and calculate diversity of this across the species present in a specific community.]

12 Issues

  • The scales of the environmental change factors are different, and so some standardisation will be needed.

13 Left-overs

Also imagine that we are studying an ecosystem that can experience simultaneous change in temperature and salinity. For example, at a particular time, the temperature might be increasing by 0.1 C per day and salinity increasing by 0.2 parts per thousand per day.

Here some text that Sam wrote for a previous ms, that we removed from there, though please confirm this before using it in here The framework proposed herein allows derivation of environment-dependent performance responses. Yet, when species are exposed to multiple environmental stressors simultaneously—a situation common in empirical systems (Bowler et al. 2020)—it is not trivial to derive performance responses. Conceptually, the principles of measuring response diversity should still apply in multiple dimensions, but the analytical formulation must necessarily differ. If multiple environmental stressors act additively, then it should be mathematically straightforward to produce a multivariate response surface with performance varying as a function of two (or more) environmental conditions. However, where multiple stressors interact, nonadditivity adds additional layers of complication to the formulation of multivariate response surfaces (e.g., Yang et al. 2022). In such cases, a different analytical approach to estimating multivariate response diversity is needed. Future work should extend the univariate response diversity principles we suggest here into a multivariate framework for measuring response diversity to multiple environmental stressors. In doing so, multivariate response diversity may link conceptually to the study of co-tolerance to multiple stressors and stress-induced community sensitivity (see Vinebrooke et al. 2004) and should prove more operationalizable than univariate approaches given the propensity for multiple environmental stressors to not only co-occur but to interact nonadditively in nature (Dieleman et al. 2012).